Statistical kinetic treatment of relativistic binary collisions 
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Abstract 

In particle-based algorithms, the effect of binary collisions is commonly described in a statistical 
way, using Monte Carlo techniques. It is shown that, in the relativistic regime, stringent constraints 
should be considered on the sampling of particle pairs for collision, which are critical to ensure 
physically meaningful results, and that nonrelativistic sampling criteria (e.g., uniform random 
pairing) yield qualitatively wrong results, including equilibrium distributions that differ from the 
theoretical Jiittner distribution. A general procedure for relativistically consistent algorithms is 
provided, and verified with three-dimensional Monte Carlo simulations, thus opening the way to 
the numerical exploration of the statistical properties of collisional relativistic systems. 

PACS numbers: 02.70.Uu, 03.30.+p, 52.65.Pp, 52.65.Rr 
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The computer-assisted kinetic analysis of the behavior of many-particle systems is funda- 
mental in several areas of modern physics, ranging from astrophysics (evolution of cosmolog- 
ical systems, dark-matter dynamics) |l| to the physics of space and laboratory plasmas (rel- 
ativistic shocks, spacecraft shielding, laser- and plasma-based particle acceleration, inertial 
confinement fusion) When the effect of close encounters (collisions) can be neglected, 
and particles can be assumed to interact via smoothly varying, long-range forces, kinetic 
particle-mesh algorithms are effective and versatile tools to study the evolution of the 
phase-space distribution function of each particle species in the system. Important examples 
are the particle-in-cell (PIC) method , which provides a self-consistent description 

of the kinetics of collisionless plasmas over distances much longer than the Debye screening 
length (as described by the Vlasov- Maxwell set of equations), and hybrid methods, mixing 
kinetic and fluid approaches . However, in situations where the inclusion of collisional 
effects in the model is critical, or when dealing with collisional-dominated many-body sys- 
tems, the collision processes must be dealt with using physically consistent algorithms, in 
order to provide a correct description of the relevant statistical properties. 

A common way to include the effect of binary collisions in particle-based algorithms (cf. 



Ref. 



m 



) is by locally changing the momenta of a suitable statistical sample of particle 



pairs using Monte Carlo (MC) techniques j9|. This approach, often referred to as Direct 
Simulation Monte Carlo (DSMC) method, provides an accurate solution of the Bolzmann 

0, 0, [n 



equation 



12l | (which is va 



ployed in molecular gas dynamics [9|, 



id for dilute systems), and has been success: 
and plasma physics 0, Q, 
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19|, 



mostly in the nonrelativistic regime. The application of these techniques to situations where 
the particle velocities are relativistic is relevant to many scenarios in High-Energy-Density 
Science, such as fast ignition of fusion targets (cf. Ref. 20|), fast electron transport in 
solid targets, proton acceleration, or shocks. As discussed in this Letter, the extension to 
relativistic regimes can not be achieved merely by guaranteeing energy-momentum conser- 
vation. Indeed, special relativity imposes further constraints on the way particle pairs are 
chosen for collision, independently of the particular type of collision process considered, 
even when the microscopic dynamics of each collision is modeled correctly. Overlooking 
these constraints on pair selection leads to unphysical results, with consequences as extreme 
as the systematic appearance of qualitatively wrong equilibrium distribution functions and 



energy-temperature relations (cf. Refs. 2l|, |22] and discussions in Refs. 23l. l24l|). 
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In this Letter, the general procedure for the statistical kinetic treatment of binary colli- 
sions in the relativistic regime is described, thus providing a consistent framework for the 
exploration of relativistic many-particle systems with MC simulations. Results from three- 
dimensional (3D) ultrarelativistic MC simulations with > 10^ computational particles are 
presented, illustrating the technique and reproducing the correct equilibrium distribution 
function over several orders of magnitudes in en ergy and particle number (Fig. [1]). A sys- 



tematic origin of conflicting results 
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2J] is identified, and, within the present 



kinetic framework, a sirnple interpretation of the numerical results published in a recent 
Letter by Cubero et al. 23| is given. 



In relativistic kinetic theory 2J, |25|, l26|| , the number of collisions, AA^, occurring within 
the space-time element AxAt about (x, t) between particles a, having momenta in the range 
(Pa, Pa + Apa), and particles 6, having momenta in the range (pb, ph -|- Ap;,), is 



AA^=^(va,Vfc)/a(x,Pa,t)/b(x,pb,i(;)Ap„Ap6AxAt, 



(1) 



where fa and are the distribution functions of species a and h (assumed to be smooth 



enough to neglect differences in the space-time coordinates before and after collisions |25| ) . 
and where A (v^, v;,) determines the collision probability as a function of the velocities and 
V;,. Since AA^ is a relativistic invariant, and so are fa, fin AxAt, then A (v^, v;,) Ap^Apb, 



and hence A (v^, v;,) •jalb, must be invariant as well [27|, with 7^ 



a,b) 



^1/2 



[a. system 



of units where the speed of light is unitary is adopted). Introducing the total cross section 
a{vy) yields A (v^, v^) '^alb = v^a{v^){l — v'^)~^/'^, which leads to the general expression 



•A. (v„, Vfc) = v,a{v,) (1 - v„ ■ Vfc) 



(2) 



where = [(v^ — v;,)^ — (v^ x V5)^]^/^/(l — ■ v^) is the absolute value of the relative 
velocity in a reference frame where one particle is at rest 
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One important consequence of Eq. ([2]) is that A(va,^b) cannot be a function of 
a single invariant parameter (e.g., v^), as it is in the nonrelativistic regime, where 
-^NR(va,Vb) = cr(|va — Vf,|)|va — Vf,|. Any violation of this essential constraint in cal- 
culations or simulations breaks the invariance of AA^, leading to qualitatively unphysical 
results, notably to equilibrium distribution functions differing from the stationary solu- 
tion of the Boltzmann equation, which, for relativistic systems, is the Jiittner function 



fj (x, p) oc exp{ru[U • p - eo7(p)]//cBT} 
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28|, 



29|, where eo is the rest energy, the 



constant U is the equilibrium mean velocity 25|, Fu = — U^, ks is the Boltzmann con- 



stant, and the invariant constant T is the equilibrium temperature measured in the reference 
frame where U = 0. 

In the statistical treatment of relativistic binary collisions, it is mandatory to adopt a 
procedure that satisfies not only the fundamental conservation laws, such as the conservation 
of the total 4-momentum, but also the relativistic invariance of AA^, a subtler but equally 
important requirement. In the nonrelativistic regime, this is usually not a concern, because 
the Galilean invariance of AA^ is trivially satisfied, with all quantities in Eq. ([T]) being 
invariant. Thus, particular attention is needed whenever applying nonrelativistic approxi- 
mations, since these may violate the invariance of AA^: a striking, paradigmatic example 
is the assumption of a uniform collision probability, ^(va,v;,) = Constant, corresponding 
to random pairing in MC algorithms, as commonly employed in nonrelativistic or weakly 
relativistic PIC simulations [7,[l4,[2^. According to special relativity, such an assumption 
is unphysical, leading to a wrong equilibrium distribution, described by a modified Jiittner 
function, /^j (x, p) oc exp{ru[U-p-eo7(p)]//i;BrMj}/{ru[eo7(p)-U-p]}, and, contextually, 
to a wrong equilibrium temperature, T^j. In the recent literature, /mj has been proposed 



as a plausib 



e extension of the nonrelativistic Maxwell-Boltzmann distribution to relativis- 



tic systems |2l|, l22|], but this possibility has been recently ruled out using one-dimensional 



(ID) numerical simulations 23|]. As shown here, /mj is obtained in MC algorithms when- 
ever ^(va,v;,) is erroneously assumed to be a function of the single invariant parameter 
fr, e.g. with uniform random pairing, independently of the particular choice of a{v,-) [as in 
the nonrelativistic case, o"(fi.) merely affects relaxation processes, having no effect on the 
equilibrium distribution]. 

In order to obtain the correct physical results in particle-based kinetic algorithms, with 
an MC approach, it is sufficient to adopt a three-step procedure: given a collection of 
particles contained in a spatial region Ax, the momenta of a statistical sample of particle 
pairs undergoing a given collision process are updated over a time interval At by 

1. Sampling the colliding pairs according to the relativistic expression of A (v^, v;,) given 



in Eq. ([2]) (e.g., with standard rejection methods [9|]), so as to guarantee the invariance 
of AA^. Depending on the problem, the number of colliding pairs must be chosen 
appropriately, ensuring that, on average, the correct number of collisions is performed, 
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and the correct collision frequency is recovered 



2. Deciding the output of each collision, using the differential cross section [24, l25|, |26| 
to evaluate the scattering angle, so as to guarantee that the microscopic details of the 
collision process are modeled correctly. For inelastic collisions (e.g., reactions, ioniza- 
tions, recombinations, pair creation/annihilation), this may involve particle generation 
and removal. 

3. Updating the momenta of all particles resulting from each collision, obeying to the 
relevant conservation laws (e.g., the conservation of the total energy-momentum and 
of the total electric charge). As an example, for an elastic collision between particles 
a and 6, this step is conveniently performed by transforming and p;, to the center- 
of-momentum frame, rotating the momenta by the appropriate scattering angle, and 



transforming the new momenta back to the laboratory frame 20 1. 



This procedure provides a correct description of the coUisional dynamics, as predicted by 
the Boltzmann equation [l^, and correctly yields the equilibrium distribution function /j, 
independently of the specific cross section, and for all energy ranges. 

As a test for the algorithm, the evolution to equilibrium of many-particle systems in 
conditions ranging from nonrelativistic to ultrarelativistic regimes has been investigated with 
massively parallel, 3D MC simulations based on the Osiris 2.0 framework 2|, employing up 
to 10^ computational particles. In the example shown here, the equilibrium distribution of a 
single species of ultrarelativistic particles undergoing elastic, isotropic collisions is analyzed 
using 2 X 10^ computational particles. A monoenergetic initial distribution has been used, 
/(x, p, t = 0) oc 5[7(p) — 7o] with 70 = 10"^ and mean velocity U = 0, and a{y^) (xl/v-s- has 
been assumed, thus yielding A (v^, v^,) oc (1 — v„ ■ v^). 

In order to provide a clear evidence that the equilibrium distribution /gq (x, p) accurately 
reproduces /j (x, p), the corresponding energy distribution 

Peq(7) = jj^ {\l ^0 + - eo7^ /eq (x, p) dpdx (3) 

has been constructed directly from the numerical data (by counting the number of particles 
having energy within finite intervals on the 7 axis), and plotted over a wide range of 7, 
spanning several orders of magnitudes (Fig. [1]). The simulated equilibrium energy distribu- 
tion accurately reproduces the theoretical curve p(7) = 71/7 — 1 exp {—e^'j/kBT), obtained 



by replacing /eq with fj in Eq. ([3]), where the equihbrium temperature fc^T = 3.33 x lO^eo 
is calculated from the initial mean energy as (7) ^ 1 + SksT/eo, the ultrarelativistic limit 
of the energy-temperature relation 



(7) = — jr- = ) (- , (4) 

eo///.dpdx K,[^) ^0 



where Kn denotes the nth order modified Bessel function of the second kind j30|. The 
numerical results are in complete agreement with the theoretical curve, correctly reproducing 
variations spanning eight orders of magnitudes in peq (Fig. [1]). 

The shape of Peqil) obtained by (incorrectly) sampling the colliding pairs according to 
the nonrelativistic approximation ^(va,v;,) = Constant, is also shown. The distribution 



reproduces the modified curve Pmj(7) = Vl ~ lexp (— eo7/^B^Mj), obtained by replacing 
/eq with /mj in Eq. (j3]). The equilibrium temperature ksT^j = 5 x lO^eo is calculated from 
the initial mean energy as (7)^^, ~ 1 + which is the ultrarelativistic limit of the 

modified energy-temperature relation 



eo + P^Mjdpdx K2 



(7)mj = '''' ' rr = T^VT' 

eo / / /Mjdpdx Ki (^^J^j j 



Although still complying with the energy- momentum conservation (step 3 above), this result 



is unphysical, because it violates t 
Lorentz-boosted reference frame 



le relativistic invariance of AA^: if performed within a 
31| . with boost factor 7b, the same simulation would exhibit 
an artificial increase of the total number of collisions by a factor 7^, as can be readily verified 
from Eqs. ([1]) and ([2]), thus yielding a significantly different dynamical evolution of the 
system and a wrong equilibrium state. 

The present analysis also allows for a straightforward kinetic interpretation of the nu- 
merical results recently presented in Ref . [23^ , where molecular dynamics (MD) simulations 
of a ID system composed of two species of colliding particles have been used to provide a 
numerical confirmation that the equilibrium one-particle distribution of a dilute relativistic 
gas is described by the Jiittner fuction /j, as opposed to the modified Jiittner function, 
/mj. The ID system considered in Ref. |23|] is a collection of impenetrable particles un- 
dergoing binary collisions, wherein interactions are zero-range and particles act as infinitely 
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FIG. 1: (Color online) Equilibrium energy spectra Peq(7) obtained using the relativistically con- 
sistent law A(ya,'Vb) oc (1 — Vq • V5) (dark) and the nonrelativistic approximation ^(va,v;,) = 
Constant (light). Markers: simulation results; solid lines: ^(7) and Pmj(7) as in the text. The 
numerical data have been taken at a single time instant after approximately 10 — 20 collisions per 
particle have occurred. Units are normalized so that J pcq{'y)d'y = 1. 

extended rigid sheets. Each collision is a localized event in space-time, with a fully deter- 
ministic outcome. Between collisions, particles are free-streaming, with the Hamiltonian 
of the system being the linear superposition of the relativistic Hamiltonians of each free 



particle, ^^(x, p) = a/cq + for the nth particle. This allows for a fully deterministic 
numerical solution of the equations of motion via standard MD techniques 321]. In the ki- 



netic approximation, the basic statistical properties of the system analyzed in Ref. 



23| (i.e., 



the equilibrium one-particle distribution function integrated over space) can be investigated 



using the relativistic Boltzmann equation, whose stationary solution is /j [2^, |25|, |26|, |29| . In 
ID, Eq. ([l]) reduces to AN = V{vr:)\va - Vb\fa{x,Pa,t) fb{x,pb,t) ApaApbAxAt, where 
V{vr) is the probability for an a-h encounter to result in a collision, with the limit 
'P(fr) 1 corresponding to impenetrable particles, as considered in Ref. 23|]. Deter- 
mining each collision event exactly, the MD algorithm used in Ref. 23(] implicitly guar- 
antees the invariance of AA^, thus yielding the correct distribution function /j. Statis- 
tical approaches recover the same result, independently of the particular shape of 'P(fr), 
provided that colliding pairs are sampled according to the relativistically consistent law 
'P{yr)\va — Vb\ = VT.V{vr){l — VaVb). As in the 3D case (Fig. [1]), if the coUiding pairs are 
erroneously sampled according to a nonrelativistic, one-parameter law of the form 'P(fi.)fr, 
the modified function /^j is always obtained. The formal proof is straightforward: in ID, the 



collision integrals [27| expressing the net change per unit time in the distribution function 
of particles a and b due to collisions read Ja,b = J 'P{vr)\va — Vb\{fafb ~ fafb)dpb,a, where 
fab ~ fa,b{x,p'^f^,t), with j)^ denoting momenta^ after collisions. Setting the local entropy 
production s{x,t) oc -Y.a=a,b I ^^sifa)Jadpa [25] to zero, then yields f'JI, = fafb, leading, 
for both species, to the equilibrium distribution function /,, with same equilibrium temper- 
ature T. If the calculation is repeated after replacing V{vr)\va — Vb\ with the nonrelativistic 
law V{v,)v„ the colhsion integrals become Ja,b = J 'P{vr)\va - VbU'j'Jalbfb - lafalbfb)dpb,a, 
yielding ■Jafalifb — lafalbfb, which leads to the modified equilibrium distribution /^j, with 
a modified equilibrium temperature Tmj. Hence, from a purely mathematical point of view, 
/mj could be considered as the stationary solution of a modified relativistic Boltzmann equa- 



tion (cf. Conclusions in Ref. 2J] and references therein), with collision integrals of the form 
Ja,b- Again, such an equation would violate the relativistic invariance of AA^, thus being 
physically inconsistent. 

In summary, the problem of providing a consistent statistical description of relativistic 
binary collisions in dilute many-particle systems has been analyzed using the standard rel- 
ativistic kinetic theory, showing that rigorous constraints hold on the way particle pairs are 
chosen for collision, and that nonrelativistic approximations (such as a uniform collision 
probability) are forbidden. By breaking the relativistic invariance of the number of collision 
events in a given space-time region, these approximations lead to unphysical, conflicting 
results, notably modified equilibrium distribution functions. Thus, in any calculation or 
simulation based on statistical sampling of colliding particles, the invariance of AA^ consti- 
tutes a fundamental validity criterion, as important as the more obvious energy-momentum 
conservation, in order to guarantee that results are physically meaningful, with the equilib- 
rium distribution being described by the Jiittner function, as predicted by the relativistic 
Boltzmann equation. The present discussion thus provides the framework for the detailed 
exploration, via Monte Carlo simulations, of the statistical properties of multi-dimensional, 
coUisional systems in the relativistic regime. 
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